function x = lsaenh(y,framesize,alpha,autoflag,ratio_or_n, Fs, start_time, end_time)
%  __  __                            _   _       
% |  \/  |                          | | | |      
% | \  / | __ _ _ __ __ _ _   _  ___| |_| |_ ___ 
% | |\/| |/ _` | '__/ _` | | | |/ _ \ __| __/ _ \
% | |  | | (_| | | | (_| | |_| |  __/ |_| ||  __/
% |_|  |_|\__,_|_|  \__, |\__,_|\___|\__|\__\___|
%                      | |                       
%                      |_|                       
%  _    _       _                    _ _         
% | |  | |     (_)                  (_) |        
% | |  | |_ __  ___   _____ _ __ ___ _| |_ _   _ 
% | |  | | '_ \| \ \ / / _ \ '__/ __| | __| | | |
% | |__| | | | | |\ V /  __/ |  \__ \ | |_| |_| |
%  \____/|_| |_|_| \_/ \___|_|  |___/_|\__|\__, |
%                                           __/ |
%                                          |___/ 
% Last Modified by Dexter Delfrate for 
% Marquette University's Speech and Signal Processing Lab.
%
% This file is part of the Marquette University filters add-on for XBAT.
%
% To run this program, open xbat.m and select 'Ephraim Malah' from filters/Marquette.
% 
%
% Ephraim Malah log spectral amplitude estimation
% based off original function "yarmal" by Ephraim & Malah
% This version includes an option to automatically find the lowest
% energy frames and use them as silence estimate.
% Author: MTJ 2009/06/29
%
% Inputs:
%       y= noisy signal
%       framesize = framesize in points (default 256)
%       alpha = smoothing constant (default 0.98)
%       autoflag = flag for auto silence-detect, 1=use (default 1)
%       ratio_or_n = if autoflag=1, silence ratio (default 0.05)
%                    if autoflag=0, noise waveform to use (No default)
%       Fs = sampling rate of the signal
%       start_time, end_time = the selection event specifying the noise waveform to use
%              
% Outputs:
%       x = enhanced signal
%
%   See also
% 
%   Copyright 2009 Marquette University, Speech and Signal Processing
%   Laboratory.



if (nargin <2)
    framesize=256;
end
if (nargin <3)
    alpha=0.98;
end
if (nargin <4)
    autoflag=1;
end
if ((nargin <5) && (autoflag==1))
    ratio_or_n=0.05;
end
if (autoflag==1)
    silratio=ratio_or_n;
elseif (autoflag==0)
    noisesig=ratio_or_n;
else
    error('Autoflag must be zero or one.');
end

y = y(:)';               % ROW vector
duration = end_time - start_time;
% framing
win_len   = framesize;
step_size = floor(framesize/2);
orig_len = length(y);
num_win = ceil((orig_len-(win_len-step_size))/step_size);
data_len = num_win*step_size+(win_len-step_size);
y=[y zeros(1,(data_len-orig_len))];  % Append extra zeros to fill out frames
x=zeros(1,data_len);
frame_start = ceil(start_time*Fs/win_len +1);
% lambdaD Noise spectrum calculation
if (autoflag==1)
    % Automatic silence detection option
    % Note silence frames are NOT overlapped - not matching signal framing
    silfr=findsil(y,win_len,silratio); % bottom 5% of energy framess
    num_win_noise=length(silfr);
    for I=1:num_win_noise
        An(:,I)=y((1+(silfr(I)-1)*win_len):(silfr(I)*win_len))';
    end
    N = fft(An.*(hanning(win_len)*ones(1,num_win_noise)));
    lambdaD = mean((N.*conj(N))');
else
    % Noise waveform option; again frames NOT overlapped
    num_win_noise=floor(duration*Fs/win_len +1);
    for I=1:num_win_noise
        An(:,I)=y((frame_start+(I-1)*win_len):frame_start -1 + (I*win_len))';
    end
    N = fft(An.*(hanning(win_len)*ones(1,num_win_noise)));
    lambdaD = mean((N.*conj(N))');    
end

% Main Loop
y = y';
for I=1:num_win
   Ay = y(1+(I-1)*step_size : win_len + (I-1)*step_size);
   Y  = fft(Ay.*hanning(win_len));
   gamaK = Y.*conj(Y)./lambdaD';
   if I==1,
      % Initial estimation of apriori SNR
      ksi=alpha+(1-alpha)*p(gamaK-1);
      TEMP=ksi./(ksi+1);
      Agal=TEMP.*exp(0.5*ei(TEMP.*gamaK)).*abs(Y);
   else
      % Estimation of apriori SNR
      ksi=alpha*Agal.*Agal./lambdaD'+(1-alpha)*p(gamaK-1);  
      TEMP=ksi./(ksi+1);
      % Estimation of the log-spectral amp. estimator
      Agal=TEMP.*exp(0.5*ei(TEMP.*gamaK)).*abs(Y);
   end
   % Adding the noisy phase
   ygal=Agal.*exp(i*angle(Y));
   % Transformation to Time-Domain; ovlap-add resynthesis
   xtag=ifft(ygal);
   x((1+(I-1)*step_size):(win_len+(I-1)*step_size))=...
   x((1+(I-1)*step_size):(win_len+(I-1)*step_size))+...
   xtag';
end
x = real(x);
x = x(1:orig_len);  % Truncate to original length signal

return


function frmind = findsil(x,frsize,p)
% function frmind = findsil(x,frsize,p)
% quick function to identify most likely silence regions in a waveform
% Note this is NOT an endpoint detector
% It just finds the lowest energy frames and concatenates them
% Note it could be improved a lot, to have minimum duration sections,
% to look at spectral content/flatness, etc.
%
% Inputs: x - signal
%         frsize - framesize in points
%         p - prior prob of noise (% of signal to use)
% Outputs: sil - concatenated silence frames
%          frmind - index of frames used 

numfr=floor(length(x)/frsize);
numsil=max(1,floor(numfr*p));

for i=1:(numfr-1)   % last one might be zero padded, don't includes
    wind=(1+(i-1)*frsize):(i*frsize);
    en(i)=sum(x(wind).^2);
end

[sen,enind]=sort(en);
frmind=sort(enind(1:numsil));  % put contiguous frames adjacent

function f=p(x)
% function f=p(x)
% Positive or not (p(x)={x,if x>0; 0,if x<0}
f=(x>0).*x;

function y=ei(x)
% ei - is an exponential integral
% ================================
% Usage:
%       y=ei(x)
% 
% Defined by :
%                inf
%  E1(x)=-Ei(-x)= S exp(-t)/t dt
%                 x

%This code caused problems during implementation. The below call to the
%   built-in exponential integral function has proven to be MUCH faster as
%   well.
% % if min(x)>0,
% %     x=(x<20).*x+(x>=20).*20;
% %     az=1;
% %     s=0;
% %     for i=1:100,
% %         az=az*i;
% %         s=s+((-x).^i)/(i*az);
% %     end;
% %     s=-(0.577215+log(x)+s);
% %     y=sign(20-x).*(sign(20-x)+1).*s/2;
% % else
% %     errordlg('the_value_cannot_be_zero_or_negative');
% % end
y = expint(x);


